######################################################################################################################
#Figure A_1: Blockwise treatment assignment
#This figure shows the assignment of districts within Jharkhand to study (10) and non-study (14) status, and the assignment of blocks within these districts to treatment and control
######################################################################################################################


################################
# Color Set up
################################

cColor <- "#97B8DE"
tColor <- "#4669BD"
eColor <- "white"
nColor <- "grey60"


################################
# Randomization map
################################


#Block shapefile
gdat.sub <- rgdal::readOGR(dsn=path.expand(paste0(AdminDataDir, "Map")), layer="SUBDISTRICT_11")

gdat.sub <- gdat.sub[gdat.sub@data$STATE_UT == "Jharkhand",]

temp <- gdat.sub@polygons[[1]]@Polygons[[1]]@coords

#align district names
gdat.sub@data$DISTRICT[gdat.sub@data$DISTRICT == "Saraikella-Kharsawan"] = "Seraikela-Kharsawan"

#Treatment status
treatment <- xlsx::read.xlsx(paste0(AdminDataDir,"Map", "/", "JH_ePOS_Blocks.xlsx"), sheetIndex = 1, stringsAsFactors = FALSE)
treatment <- treatment[1:132,c("District","Block","Treated")] #remove extra row & uneccessary columns

#align district names to those in shapefile
treatment$District[treatment$District == "Saraikella-Kharsawan"] <- "Seraikela-Kharsawan"
treatment$DISTRICT <- treatment$District
treatment$District <- NULL
treatment$NAME <- treatment$Block
treatment$Block <- NULL


#Merge data
gdat.sub@data$order <- row.names(gdat.sub)
#must use join here; if use "merge", then data will be scrambled and matched to incorrect polygons
gdat.sub@data <- plyr::join(gdat.sub@data,treatment, by = c("DISTRICT","NAME"),type = "left", match = "all") 
gdat.sub@data$Treated[is.na(gdat.sub@data$Treated)] <- -1


#create a block data frame for ggplot
gdat.sub.points <- fortify(gdat.sub, region = "C_CODE11")
gdat.sub.points <- merge(gdat.sub.points,gdat.sub@data, by.x = "id", by.y = "C_CODE11")


#create a district dataframe for ggplot
gdat.dist <- maptools::unionSpatialPolygons(gdat.sub,gdat.sub@data$DISTRICT)
gdat.dist.points <- fortify(gdat.dist, region = "DISTRICT")


################################################################
# Output Plot
################################################################

ggplot(data = gdat.sub.points) + 
  geom_polygon(aes(x = long, y = lat, group = id, fill = factor(Treated))) + 
  geom_path(aes(x = long, y = lat, group = id), color = "grey60") + 
  geom_path(data = gdat.dist.points ,aes(x = long, y = lat, group = group), color = "grey30") + 
  scale_fill_manual(values = c(eColor, cColor, tColor), labels = c("Non-study  ", "Control  ", "Treatment"))+
  theme(panel.background = element_blank(), legend.position = c(.8, .05),
        legend.title = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_blank(), 
        axis.text.x = element_blank(), axis.text.y = element_blank(),
        axis.ticks = element_blank(), axis.title.x = element_blank(), 
        axis.title.y = element_blank(), legend.text=element_text(size=15),
        legend.key = element_rect(fill = "black", color = NA),
        legend.background = element_blank(),legend.direction="horizontal",
        legend.box.background = element_rect(colour="white"))

while (!is.null(dev.list()))  dev.off()

ggsave(file=paste(OutputDir, "FigureA_1.pdf", sep=""), width=12, height=12)

